1 Overview

This document contains the code from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code includes the case study analysis on ringed seal movement in eastern Hudson Bay, Canada.

2 Prepare

2.1 Load packages

library(tidyverse)
library(ggplot2)
library(lubridate)
library(rnaturalearth)
library(rgdal)
library(terra)
library(amt)
library(momentuHMM)
library(viridis)
library(sf)
library(here)

2.2 Fish data

Next we load in the fish data. This is a subset of the data from Florko et al. (2021).

# load fish data
fish <- read.csv(here("data/preydiv.csv"))
head(fish)
##        lon     lat preydiv
## 1 168329.9 -116109   0.513
## 2 168329.9 -166109   0.503
## 3 168329.9 -216109   0.512
## 4 168329.9 -266109   0.533
## 5 168329.9 -316109   0.557
## 6 168329.9 -366109   0.575

Visualize the fish data.

#prepare world data
## List of azimuthal equal area - HudsonBay 
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string
## Warning in wkt(obj): CRS object has no comment
# plot fish map
fishmap <- ggplot() + 
  geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
  scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 
fishmap

2.3 Seal data

Next we load in the seal movement data. This is a subset of the data from Florko et al. (2023).

seal <- read.csv(here("data/seal_track_m.csv")) 
head(seal)
##        lon       lat             date       id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1
# ensure the data is in the correct format
seal <- seal %>%
  mutate(id = as.character(id),
         date = as.Date(date)) 

Visualize the seal data on top of the fish data.

# plot
sealfishmap <- fishmap + 
  geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") + 
  geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE") 

sealfishmap

3 Fit models

3.1 RSF

Fit using amt (Signer et al. 2019)

# generate availability sample
set.seed(2023)
data_rsf <- seal %>%
  make_track(lon, lat, date) %>%
  random_points() # default is ten times as many available points as observed points

# plot used vs available locations
sealfishmap +
  geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
  scale_color_manual(values = c("black", "#FCEEAE"), 
                     label = c("Available", "Observed"), name = "Data type") 

# rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

data_rsf <- data_rsf %>%
  extract_covariates(fish_raster)

# fit rsf
rsf1 <- data_rsf %>%
  amt::fit_rsf(case_ ~ preydiv, model = TRUE)

# see parameter estimates
summary(rsf1)
## 
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"), 
##     data = data, model = TRUE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.5753  -0.4729  -0.4294  -0.3754   2.3363  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -6.118      1.400  -4.369 1.25e-05 ***
## preydiv        6.353      2.312   2.748    0.006 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 938.28  on 1539  degrees of freedom
## Residual deviance: 930.60  on 1538  degrees of freedom
## AIC: 934.6
## 
## Number of Fisher Scoring iterations: 5

3.2 SSF

Fit using amt (Signer et al. 2019)

# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
  make_track(lon, lat, date) %>%
  steps() %>%
  random_steps() %>% # default is ten times as many available points as observed points
  arrange(case_)

# plot used vs available locations
sealfishmap +
  geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
  scale_color_manual(values = c("black", "#FCEEAE"), 
                     label = c("Available", "Observed"), name = "Data type") 

# extract prey diversity covariate
data_ssf <- data_ssf %>%
  extract_covariates(fish_raster, where = "end") #sample at end of step

# transform movement covariates
data_ssf <- data_ssf %>%
  mutate(cos_ta = cos(ta_),
         log_sl = log(sl_))

# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
  fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)

## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
  fit_clogit(case_ ~ log_sl * preydiv + cos_ta + strata(step_id_), model = TRUE)

# see parameter estimates
summary(ssf1)
## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl + 
##     cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
## 
##   n= 1516, number of events= 138 
##    (2 observations deleted due to missingness)
## 
##              coef exp(coef)  se(coef)      z Pr(>|z|)
## preydiv  0.957799  2.605953  6.772711  0.141    0.888
## log_sl  -0.008404  0.991631  0.083328 -0.101    0.920
## cos_ta   0.031707  1.032215  0.130643  0.243    0.808
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## preydiv    2.6060     0.3837 4.477e-06 1.517e+06
## log_sl     0.9916     1.0084 8.422e-01 1.168e+00
## cos_ta     1.0322     0.9688 7.990e-01 1.333e+00
## 
## Concordance= 0.494  (se = 0.029 )
## Likelihood ratio test= 0.09  on 3 df,   p=1
## Wald test            = 0.09  on 3 df,   p=1
## Score (logrank) test = 0.09  on 3 df,   p=1
summary(ssf2)
## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ log_sl * preydiv + 
##     cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
## 
##   n= 1516, number of events= 138 
##    (2 observations deleted due to missingness)
## 
##                      coef  exp(coef)   se(coef)      z Pr(>|z|)
## log_sl          1.663e+00  5.273e+00  1.286e+00  1.293    0.196
## preydiv         2.791e+01  1.318e+12  2.177e+01  1.282    0.200
## cos_ta          2.338e-02  1.024e+00  1.311e-01  0.178    0.858
## log_sl:preydiv -2.744e+00  6.433e-02  2.101e+00 -1.306    0.192
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## log_sl         5.273e+00  1.897e-01 4.241e-01 6.556e+01
## preydiv        1.318e+12  7.589e-13 3.919e-07 4.430e+30
## cos_ta         1.024e+00  9.769e-01 7.917e-01 1.324e+00
## log_sl:preydiv 6.433e-02  1.555e+01 1.046e-03 3.955e+00
## 
## Concordance= 0.555  (se = 0.029 )
## Likelihood ratio test= 1.81  on 4 df,   p=0.8
## Wald test            = 1.81  on 4 df,   p=0.8
## Score (logrank) test = 1.81  on 4 df,   p=0.8

3.3 HMM

Fit using momentuHMM (McClintock and Michelot, 2018)

# prepare dataset
data_hmm <- seal %>%
  make_track(lon, lat, date) %>%
  extract_covariates(fish_raster) %>%
  mutate(ID = 1,
         x = x_,
         y = y_,
         date = t_) %>%
  dplyr::select(ID, x, y, date, preydiv) %>%
  as.data.frame()

# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"), 
                                 type = "UTM",
                                 covNames = c("preydiv"))

# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution

mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000,  10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0) 
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state

formula = ~ preydiv # identify covariates

# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
            dist=list(step=stepDist,angle=angleDist),
            Par0=list(step=stepPar0,angle=kappa0))

# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)

# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
            dist=list(step=stepDist,angle=angleDist),
            Par0=Par0_hmm1$Par,
            delta0 = Par0_hmm1$delta, 
            beta0 = Par0_hmm1$beta,
            formula=formula)

# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() + 
  scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
  geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) + 
  geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) + 
  scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), 
                     labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"), 
                     name = "HMM state") +
  scale_shape_manual(values = c(15, 16, 17), 
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"), 
                     name = "HMM state") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plot

#---------------------PLOT HISTOGRAMS

## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
  make_track(lon, lat, date) %>%
  extract_covariates(fish_raster) %>%
    mutate(ID = 1,
         x = x_,
         y = y_,
         date = t_) %>%
  dplyr::select(ID, x, y, date, preydiv) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("x", "y")) %>%
  st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
  st_transform("+proj=longlat +datum=WGS84") %>%
  mutate(long = unlist(map(geometry,1)),
          lat = unlist(map(geometry,2))) %>%
  st_drop_geometry()

# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))


# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution

mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)

# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
            dist=list(step=stepDist,angle=angleDist),
            Par0=list(step=stepPar0,angle=kappa0))

# get parameters 
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)

# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
            dist=list(step=stepDist,angle=angleDist),
            Par0=Par0_hmm1_km$Par,
            delta0 = Par0_hmm1_km$delta, 
            beta0 = Par0_hmm1_km$beta,
            formula=formula)

# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))


##-----STEP DIST histogram

#calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)

#plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")

#generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)

#get converged mean and sd for each state 
meanARS <- hmm2_km$mle$step[1,1]  
sdARS <- hmm2_km$mle$step[2,1]    

meanCR <- hmm2_km$mle$step[1,2]    
sdCR <- hmm2_km$mle$step[2,2]     

meanTR <- hmm2_km$mle$step[1,3]    
sdTR <- hmm2_km$mle$step[2,3]    

#calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}

#get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS),  scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR),  scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR),  scale=sc(meanTR,sdTR)) * stateFreq[[3]]


#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR,  State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR,  State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)

#plot distributions
hmmstepdist_plot <- ggplot() +
  geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
  scale_colour_manual(values=c(colours.states,"#000000"), 
              breaks = c('Foraging', 'ARS', 'Travelling'), 
              labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  scale_linetype_manual(values=c("solid","solid", "solid"), 
              breaks = c('Foraging', 'ARS', 'Travelling'), 
              labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  ylab("Density") + 
  xlab("Step length (kms)") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.

## Warning: Please use `linewidth` instead.
hmmstepdist_plot

##-----ANGLE DIST
# #generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)

#get converged mean and concentration for each state 
meanARS <- hmm2_km$mle$angle[1,1]  
sdARS <- hmm2_km$mle$angle[2,1]    

meanCR <- hmm2_km$mle$angle[1,2]    
sdCR <- hmm2_km$mle$angle[2,2]  

meanTR <- hmm2_km$mle$angle[1,3]    
sdTR <- hmm2_km$mle$angle[2,3]  

#get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS,  kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR,  kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR,  kappa=sdTR) * stateFreq[[3]]


#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR,  State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR,  State="Travelling",x=x)

cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)

#plot distributions
hmmangledist_plot <- ggplot() +
  geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
  scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
  scale_x_continuous(limits=c(-pi,pi))+
  ylab("Density") + 
  xlab("Turn angle (radians) ") +
  theme_minimal()
hmmangledist_plot

4 Plot prediction maps

Prep the data

# prep the fish data
newfish <- fish_raster %>%
  terra::as.data.frame(xy = TRUE) %>%
  filter(x > 100000 & x < 600000 & y > -550000 & y < 0)

4.1 RSF

modcoef <- summary(rsf1)$coef
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))
x <- x / (1 + x)

# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
newfish$rsf_prediction <-range01(x)

# plot
map_RSF <- ggplot() + 
  geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
  scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 
## Regions defined for each Polygons
map_RSF

4.2 SSF

# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
  mutate(date = as.POSIXct(date)) %>%
  make_track(lon, lat, date) %>%
  steps() %>%
  random_steps() %>% 
  arrange(case_) %>%
  amt::extract_covariates(fish_raster, where = "both") %>%
  na.omit() 

### -------------- SSF1
# fit SSF1 model
m1 <- data_ssf |> 
  fit_clogit(case_ ~ preydiv_end +  cos(ta_) + log(sl_) + 
              # add terms for home range
               x2_ + y2_ + I(x2_^2 + y2_^2) +
               strata(step_id_))


# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
                       mutate(date = as.POSIXct(date)) %>%
                       make_track(lon, lat, date))
                          # random sample along the track as the start
                          [sample(1:nrow(seal), 1), ], 
                    dt_ = hours(2)) 

# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
                            stochastic = TRUE,
                            tolerance.outside = 0.1,
                            n.control = 1e3)

# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points. 
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)

# plot simulated track
ssf_track_1 <- fishmap +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  geom_point(data = p1, aes(x = x_,y = y_), alpha = 0.61) +
  geom_path(data = p1, aes(x = x_,y = y_)) 
## Regions defined for each Polygons
# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1, 
    x_ = p1$x_, y_ = p1$y_,
    t_ = p1$t_, dt = p1$dt) |> 
  filter(!is.na(x_)) |> 
  make_track(x_, y_) |> 
  hr_kde(trast = fish_raster, which_min = "global") %>%
  hr_ud() %>% 
  terra::as.data.frame(xy = TRUE)

# plot SSUD
map_SSF1 <- ggplot() + 
  geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
  scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 
## Regions defined for each Polygons
map_SSF1

### -------------- SSF2
# fit SSF2 model
m2 <- data_ssf |> 
  fit_clogit(case_ ~ preydiv_end +  preydiv_end:cos(ta_)+
               preydiv_end:log(sl_) + 
               # add terms for home range
               x2_ + y2_ + I(x2_^2 + y2_^2) +
               strata(step_id_))


# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
                       mutate(date = as.POSIXct(date)) %>%
                       make_track(lon, lat, date))
                          # first location of the track as the start
                          [sample(1:nrow(seal), 1), ], 
                    dt_ = hours(2)) 

# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
                            stochastic = TRUE, 
                            tolerance.outside = 0.3,
                            n.control = 1e3)

# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points. 
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)

# plot simulated track
ssf_track_2 <- fishmap +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
  geom_path(data = p2, aes(x = x_,y = y_)) 
## Regions defined for each Polygons
# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1, 
    x_ = p1$x_, y_ = p1$y_,
    t_ = p1$t_, dt = p1$dt) |> 
  filter(!is.na(x_)) |> 
  make_track(x_, y_) |> 
  hr_kde(trast = fish_raster, which_min = "local") %>%
  hr_ud() %>% 
  terra::as.data.frame(xy = TRUE)

# plot SSUD
map_SSF2 <- ggplot() + 
  geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
  scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) 
## Regions defined for each Polygons
map_SSF2

library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
plot_grid(map_SSF1, map_SSF2, ssf_track_1, ssf_track_2, nrow =2)

4.3 HMM

x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3 

newfish_long <- newfish %>%
  tidyr::pivot_longer(cols = hmm_state1:hmm_state3, 
               names_to = "model", values_to = "prediction") %>%
  mutate(dplyr::across(model, factor, levels=
               c("hmm_state1", "hmm_state2", "hmm_state3")))

map_HMM <- ggplot() + 
  geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
  scale_fill_viridis(option = "mako", limits = c(0,1)) +
  labs(fill = 'HMM predicted\nprobability') +
  geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
  coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
  facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
                                              'hmm_state2' = "Area-restricted search",
                                              'hmm_state3' = "Travel")))
## Regions defined for each Polygons
map_HMM

5 Plot prey div vs probability of selection

5.1 RSF

base <- newfish %>% 
  #  (note, it doesn't matter what values you pick for sl_ and ta_)
  mutate(log_sl = log(45),
         cos_ta = cos(1))

x1 <- base 

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Report status
  cat(i, "of", nrow(x1), "\n")
  # Calculate log-RSS for that row
  xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})
## 1 of 72 
## 2 of 72 
## 3 of 72 
## 4 of 72 
## 5 of 72 
## 6 of 72 
## 7 of 72 
## 8 of 72 
## 9 of 72 
## 10 of 72 
## 11 of 72 
## 12 of 72 
## 13 of 72 
## 14 of 72 
## 15 of 72 
## 16 of 72 
## 17 of 72 
## 18 of 72 
## 19 of 72 
## 20 of 72 
## 21 of 72 
## 22 of 72 
## 23 of 72 
## 24 of 72 
## 25 of 72 
## 26 of 72 
## 27 of 72 
## 28 of 72 
## 29 of 72 
## 30 of 72 
## 31 of 72 
## 32 of 72 
## 33 of 72 
## 34 of 72 
## 35 of 72 
## 36 of 72 
## 37 of 72 
## 38 of 72 
## 39 of 72 
## 40 of 72 
## 41 of 72 
## 42 of 72 
## 43 of 72 
## 44 of 72 
## 45 of 72 
## 46 of 72 
## 47 of 72 
## 48 of 72 
## 49 of 72 
## 50 of 72 
## 51 of 72 
## 52 of 72 
## 53 of 72 
## 54 of 72 
## 55 of 72 
## 56 of 72 
## 57 of 72 
## 58 of 72 
## 59 of 72 
## 60 of 72 
## 61 of 72 
## 62 of 72 
## 63 of 72 
## 64 of 72 
## 65 of 72 
## 66 of 72 
## 67 of 72 
## 68 of 72 
## 69 of 72 
## 70 of 72 
## 71 of 72 
## 72 of 72
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)

# PLOT
rsf_line <- ggplot(res, aes(x = preydiv_x1, y = (log_rss))) +
  geom_line(size = 1, color = "#F8D59F",linetype = 2) +
  geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()

5.2 SSF

## Prediction for SSF1

x1 <- base 

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Calculate log-RSS for that row
  xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})

# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)


# PLOT
ssf_line_1 <- ggplot() +
  geom_line(data = res, aes(x = preydiv_x1, y = log_rss), size = 1, linetype = 3) +
  geom_ribbon(data = res, aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.3) +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()


## Prediction for SSF2
nums = seal %>%
  make_track(lon, lat, date) %>%
  steps %>%
  summarize(quants = quantile(sl_, c(0.25, 0.5, 0.75))) %>%
  pull()

x1 <- base 

results <- lapply(nums, function(i) {
x1$log_sl <- i

x2 <- base  %>% 
  mutate(preydiv = mean(base$preydiv))

# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.

log_rss_list <- lapply(1:nrow(x1), function(i) {
  # Calculate log-RSS for that row
  xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
  # Return the element $df
  return(xx$df)
})

# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)

} )
results <- dplyr::bind_rows(results) %>%
  mutate(log_sl_x1 = as.factor(round(log_sl_x1, 1)))

# PLOT
ggplot(results, aes(x = preydiv_x1, y = (log_rss))) +
  geom_line(size = 1, aes(color = log_sl_x1, group = log_sl_x1)) +
  geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = log_sl_x1, group = log_sl_x1), alpha = 0.3) +
  xlab("Prey diversity") +
  ylab("log-RSS") +
  theme_minimal()

# Plot
ssf_line_1 +
  geom_line(data = results, size = 1, aes(x = preydiv_x1, y = (log_rss), 
                                          color = log_sl_x1, group = log_sl_x1)) +
  geom_ribbon(data = results, aes(ymin=lwr, ymax=upr, x=preydiv_x1, 
                                  fill = log_sl_x1, group = log_sl_x1), alpha = 0.3)

5.3 HMM

# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)

state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)

# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
  mutate(state = as.character(state)) 
# plot
hmmprobs_plot <- ggplot() + 
  geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
  geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state), 
              alpha = 0.4, show.legend = TRUE) +
  ylab("Stationary state probabilties") +
  xlab("Prey diversity") +
  scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
    scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
                     labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
  theme_minimal()

hmmprobs_plot

6 References

Florko, K.R.N., Tai, T.C., Cheung, W.W.L., Sumaila, U.R., Ferguson, S.H., Yurkowski, D.J., Auger-Méthé, M. 2021. Predicting how climate change threatens the prey base of Arctic marine predators. Ecology Letters, 24: 2563-2575. https://doi.org/10.1111/ele.13866
Florko, K.R.N., Shuert, C.R., Cheung, W.W.L., Ferguson, S.H., Jonsen, I.D., Rosen, D.A.S., Sumaila, U.R., Tai, T.C., Yurkowski, D.J., Auger-Méthé, M. 2023. Linking movement and dive data to prey distribution models: new insights in foraging behavior and potential pitfalls of movement analyses. Movement Ecology, 11:17 https://doi.org/10.1186/s40462-023-00377-2
McClintock, B.T., Michelot, T. 2018. momentuHMM: R package for generalized hidden markov models of animal movement. Methods Ecol. Evolut. 9, 1518-1530. https://doi.org/10.1111/2041-210X.12995
Signer, J., Fieberg, J., & Avgar, T. (2017). Estimating utilization distributions from fitted step‐selection functions. Ecosphere, 8(4), e01771. https://doi.org/10.1002/ecs2.1771
Signer, J., J. Fieberg, and T. Avgar. 2019. Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses. Ecology and Evolution 9:880–890. https://doi.org/10.1002/ece3.4823